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Abstract. We study the problem of non-conventional Anderson localization 
emerging in bilayer periodic-on-average structures with alternating layers of materials 
with positive and negative refraction indices n a and n^. Main attention is paid 
to the model of the so-called quarter stack with perfectly matched layers (the 
same unperturbed by disorder impedances, Z a = Zb, and optical path lengths, 
n a d a = \nb\db, with d a , db being the thicknesses of basic layers). As was recently 
numerically discovered, in such structures with weak fluctuations of refractive indices 
(compositional disorder) the localization length Li oc is enormously large in comparison 
with the conventional localization occurring in the structures with positive refraction 
indices only. In this paper we develop a new approach which allows us to derive the 
expression for Li oc for weak disorder and any wave frequency ui. In the limit lo — > one 
gets a quite specific dependence, oc ct 4 u; 8 which is obtained within the fourth order 
of perturbation theory. We also analyze the interplay between two types of disorder, 
when in addition to the fluctuations of n a , rib the thicknesses d a , db slightly fluctuate 
as well (positional disorder). We show how the conventional localization recovers with 
an addition of positional disorder. 
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1. Introduction 

Due to a remarkable progress in manufacturing one-dimensional (ID) systems with 
given transport characteristics the interest to a rigorous analysis of various ID models 
has greatly increased in recent decades (see, e.g., book [1] and references therein). 
Of a particular interest are bilayer structures in optics [2] and electromagnetics [3], 
semiconductor superlattices [I], in the devices with alternating quantum wells and 
barriers in electronics, etc. 

Unlike the disordered models with continuous potentials for which the theory is 
fully developed, the analysis of the Kroing- Penney type models meets serious theoretical 
difficulties. One of the first non-trivial problems was a famous effect of band center 
anomaly in the standard tight-binding Anderson model. As was found numerically [5], 
the value of the localization length at the band center did not support simple analytical 
predictions based on the conventional perturbation theory [6]. More careful analysis 
performed in [7J |8] have shown that one has to use specific technics in order to derive 
correct results. The origin of the discovered effect was found to be related to a kind 
of resonances emerging for weak disorder. As was later understood, such an anomaly 
is due to a non-homogeneous distribution of the phase of a wave function, generated 
in the process of wave propagation along a structure. Similar effects are also known to 
occur in the vicinity of band edges of energy spectra (see, for example, [9] and references 
therein) . 

To date, it is understood that the effect of resonance transmission of a wave 
through any periodic structure poses a similar problem of a correct description of various 
transport characteristics for the models that are periodic-on-average. The latter term 
refers to various kinds of weak disorder added to a strictly periodic potential. This 
type of systems emerge naturally since in practice on the top of periodic structures such 
as photonic lattices or electron superlattices, weak random variations of parameters 
are experimentally unavoidable. Such variations can occur for the width of barriers 
or wells, electromagnetic or optical characteristics of materials, effective masses of 
electrons in superlattices, etc. For the frequency of waves far from the resonances 
the phase distribution of propagating waves is typically constant. In this case the 
analytical treatment of the transmission coefficient or localization length is relatively 
easy. However, in the vicinity of band edges or resonances the distribution of wave phase 
turns out to be highly non-homogeneous, the fact that leads to a dramatic complication 
of a rigorous analysis (see review [9]). 

Another effect which needs specific analytical tools in its analysis is the influence of 
various correlations in a disorder. As is already understood, the underlying correlations, 
either of short or long-range type, can significantly suppress or enhance the localization 
length, therefore, strongly influence the global characteristics of the transmission [9] . In 
particular, by imposing specific long-range correlations one can practically create the 
devices with given transmission/reflection characteristics. This effect of enhancement 
or suppression of the Anderson localization was recently observed experimentally in the 
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waveguides with both the bulk [TU] or surface [TT] correlated disorder. It should be 
stressed that specific long-range correlations can emerge naturally in physical systems, 
and not only in the systems with intentionally included correlations. One of such 
examples will be discussed in this paper. 

New perspectives in creating the devices with unusual transport characteristics are 
related to specific optic properties of metamaterials embedded in periodic structures 
[121 [131 HH [151 US IH riEl [TSl [2U1 EH [22]. A particular system of an increasing 
interest is an array of two alternating a- and b -layers with equal optical path lengths, 
n ad a = \nb\di,. Here the a -layers are made of right-handed (RH) material with positive 
refractive index, n a > and thickness d a . In contrast, the b -layers refer to the left- 
handed (LH) material with negative index, nj, < 0, and thickness d&. In such a model 
the phase shift of the wave gained in the a -layer is fully compensated by the subsequent 
shift in the next b -layer. As a result, the total phase shift vanishes after passing any of 
N units constructed by (a, b) cells. If, in addition, the layers are matched (Z a = Z b ), 
this results in a kind of invisibility of the structure for an observer. 

Recent numerical data [H] obtained for the array of two matched alternating RH 
and LH layers with d a = db and weakly disordered refractive indices n a m ~ 1, 
have demonstrated enormously fast divergence of the localization length, L~[ oc oc co 6 for 
u> — > 0. More detailed numerical study [15] has shown that the power k in dependence 
L^l oc u K increases with the sample length, and approaches the value k ph 8.78. This 
result is in contrast with the conventional dependence Lf ], oc u 2 known to occur for 
many models, both with continuous and periodic-on-average potentials. 

As was found in article [23], for the model discussed in [HI [15] the phase of wave 
propagating through the RH-LH array with fluctuating refractive indices is described 
by a highly non-uniform distribution. This effect is similar to that arising in the tight- 
binding Anderson model for the energy close to band edges. Further analysis [23] has 
led to a remarkable conclusion that the localization length diverges in the quadratic 
approximation in disorder strength. In the next study [23] a new method was suggested 
allowing one to resolve the above problem of non-conventional Anderson localization. 
It was shown that the inverse localization length within the fourth-order perturbation 
theory is described by the asymptotics Lj^ oc a A u s for small u. Therefore, the found 
dependencies with k«6 have to be attributed to a not large enough size N of samples 
used in numerical calculations. 

The aim of this paper is two-fold. First, we present a full description of the method 
used in [2^ when analyzing the abnormal localization, numerically observed in [TU [15] . 
We generalize the model studied in [HI [15l |23j [24] and give many details which are 
important for understanding our method, together with a comprehensive discussion of 
the obtained results. Second, we present new analytical results for the model with two 
kinds of disorder, namely, with the positional and compositional disorders, and show 
the interplay between two mechanisms of localization. We show how the abnormal 
localization reduces to the conventional one in dependence on the relative strengths of 
the two disorders. 
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2. Model formulation: Hamiltonian map 

We consider the propagation of an electromagnetic wave of frequency u in an infinite 
dielectric array (stack) of two alternating a- and b -layers. Every kind of layers 
is respectively specified by their thickness d ajb , dielectric permittivity e a , b , magnetic 
permeability fi a ,b, refractive index n a ^ = y/e a , b fi a ,b, impedance Z a ^ = fi a ,b/n atb and wave 
number k a>b = um a ^/c. We shall analyze two systems: the homogeneous stack in which 
both a- and b -layers are made of RH materials, and mixed stack where the a -layers 
contain the RH material while b -layers are composed of the left-handed (LH) material. 
As is known, in the RH medium all optic parameters are positive. On the contrary, 
in the LH medium the permittivity, permeability and, consequently, the refractive 
index are negative, however, the impedance remains to be positive. In what follows, 
we consider the situation when the disorder is originated from random variations of 
both layer thicknesses d a ^ {positional disorder), as well as from the fluctuations of e a ^ 
{compositional disorder). Specifically, we assume a weakness of both types of disorder: 

dan = d a [l + g a {n)], (d an ) = d a ; (2.1a) 

d hn = d b [l + Qb(n)], (dbn) = d b . ( 2 -16) 

and 

£ an = e a [l + Va{n)] 2 , n an = n a [l + r} a (n)], 

Z an = Z a [l + ^{n)]' 1 , k an = um a [l + rj a {n))/c; (2.2 a) 

e bn = e h [\ + r, h {n)} 2 , n hn = n b [l + rj b {n)], 

Z bn = Z b [l + ?7 6 (n)] _1 , k bn = um b [l + r) b {n)\/c. (2.26) 

Here the index n enumerates the n-th unit (a, 6) cell, and we assume that the magnetic 
permeabilities \i a and \i b are disorder-independent. 

The random sequences g a , b {n) and rj ajb {n) imposing, respectively, the positional 
and compositional disorder are specified by white-noise entries with the zero average 
and small variances a 2 . , 1 , 

(p a ,b{ n )) = °> (°l,b( n )) = v% (p a {n)g b {n')) = a 2 g 5 ab 5 nn > ; 
(Va,b{ n )) = 0> (vl, b ( n )) = rf, (Va{n)vb{n')) = o- 2 5 ab 5 nn , ; (2.3) 

(Qa,b{ n )Va,b{ n ')) = • 

Hereinafter, the angular brackets (...) stand for the averaging over different realizations 
of a random structure {ensemble averaging) or along its single realization {spatial 
averaging), which is assumed to be equivalent due to ergodicity. Numerically, when 
generating random sequences g a , b {n) and rj ajb {n) we use the flat distribution on a finite 
interval. However, our analytical expressions are valid for any distribution with finite 
variance. 
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Within every a - or 6 -layer the electric field of the wave, 

E(x, t) = E(x) exp(-iut), (2.4) 

obeys the one-dimensional Helmholtz equation with two boundary conditions at the 
interfaces between neighboring layers, 

' + kl b ) E a>b (x) = 0, (2.5) 



dx 2 

E a ( Xl ) = E b { Xi ), ^E'^x,) = fx b l E' b (x t ). (2.6) 

The x-axis is directed along the array of bilayers with x = x% standing for the interface 
coordinate. The prime implies the derivative with respect to x. 

The general solution of the wave equation (12.51) inside the n-th unit (a, 6) cell can 
be presented as a superposition of two standing waves, 

E a (x) = E a (x an )cos [k an (x - x an )] + k~^E' a (x an ) sin [k an (x - x an )] (2.7a) 

inside a n layer, where x an ^ x ^ x bn ; 
E b (x) = E b (x bn ) cos [k bn (x - x K )\ + k^E' b {x bn ) sin [k bn (x - x K )} (2.7b) 

inside b n layer, where x bn ^ x ^ x an+1 . 

The coordinates x a „ and x bn refer to the left-hand edges of successive a n and b n layers, 
respectively. The thickness of individual layers is defined as 

dan = x bn - x an and d bn = x an+1 - x bn . (2.8) 

The solutions (12. 7 aft and (12.761) give rise to useful relations between the electric field 
E a b and its derivative E' a b at opposite boundaries of a n and b n layers, 

E a (x bn ) = E a (x an ) cos if an + k^E' a (x an ) sin(p ani 

K( X b n ) = -kanE a (x an ) Sm(p an + E' a (x an ) COS '■fan] (2.9 a) 

E b (x an+1 ) = E b (x bn ) cos^ bn + k b ^E' b (x bn ) sin^ bn , 

E 'b( x a n+1 ) = -k bn E b (x bn ) sin^ bn + E b (x bn ) cos^ 6n . (2.96) 

The disordered phase shifts <p an and ip bn depend on the cell index n due to randomized 
refractive indices n an and n bn , as well as via random thicknesses d an and d bn of the 
layers, 

fan = k an d an = <p a [l + r) a (n)][l + Q a {n)], (2.10a) 
<fbn = h n d bn = <p b [l + ?7 b (n)][l + Qb{n)]. (2.106) 
Here the unperturbed phase shifts <p a and ip b are defined by the expressions: 

<-Pa = k a d a = un a d a /c, <p b = k b d b = un b d b /c. (2.11) 
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Then, the combination of relations (12. 9 aft . ( 12.961) with the boundary conditions ( 12. 6 p 
at the interfaces Xi = x bn and Xi = x an+1 , yields the recurrent relations describing the 
wave transfer through the n-th unit (a, b) cell, 

Qn+i = A n Q n + B n P n , P n +i = ~C n Q n + D n P n . (2-12) 

Here Q n and P n refer to the normalized electric field and its derivative, respectively, 
taken at left-hand edge of the n-th unit (a, b) cell, 

Q n = Z^ 2 E a (x an ), P n = (c/u)Z l J 2 E' a (x an ). (2.13) 

The normalization factors in (I2.13P contain the unperturbed impedance Z a of a -layers, 
see definitions ( 12.2 ah and ( 12.2&I) . The randomized fact ors v4 n , B n , C n -i D n read 

A n = cos ip an cos y? bn - Z m ]Z bn sin an sin (fi bn , 

B n Z a Z an Sin if an COS (f bn + Z bn COS if an sin Lf bn , 
CnZa 1 = Zan sin Van COS (p bn + COS ip an Sm (p bn , 

D n = cos (p an cos (p bn - Z an Z^ sin (p an sin ^ bn . 

As one can see, for a non-zero disorder the coefficients ( 12.141) depend on the cell index n 
due to randomized phase shifts f a n,fbn an d impedances Z an , Z bn . It should be stressed 
that the phase shifts are influenced by both the compositional and positional disorders, 
however, the impedances depend on the compositional disorder only, see relations (12. lap . 
(l2Tfej) . (E^ajl . (I2T2I and (l2TT0^) . (l2~T0fej) . Note also that the recurrent relations (l2~T2l 
belong to the class of area-preserving maps whose determinant equals one, 

A n D n + B n C n = l. (2.15) 

Remarkably, the relations ( I2.12P can be treated as the Hamiltonian map describing 
the evolution of trajectories in discrete time n for a linear oscillator subjected to a time- 
dependent parametric force. In such a representation Q n and P n can be treated as the 
classical coordinate and momentum, respectively. Therefore, the problem of quantum 
localization can be formally reduced to the analysis of the properties of trajectories in 
the phase space (Q,P). For the analytical treatment it is convenient to pass to polar 
coordinates, namely, to the radius R n and angle 9 n , 

Q n = R n cos9 n , P n = R n sm9 n . (2.16) 

It can be shown that the Hamiltonian map in the radius- angle presentation gets the 
form, 

-Rn+i\ 2 _ dfl n+ i arctan ( — ^ >n ^ an ^ n \ (2 17) 

Then the localization length L; oc can be derived as the inverse of the Lyapunov exponent 
A in accordance with the definition [9j 123] , 

d -A=i/lnf^lV\ = -i/ln^l>. (2.18 



Lioc 2 \ V Rn ) I 2 \ d# 
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Here the average (. . .) is performed along the "time" n, with a consequent average 
over different realizations of disorder in order to reduce the fluctuations. Such a 
representation is very useful since instead of the two-dimensional map for R n and 9 n the 
analysis of the Lyapunov exponent can be done in terms of the one-dimensional map for 
the phase 6 n . It should be, however, noted that the expression (12.18j) is not valid in the 
vicinity of band edges, where special methods have to be used (see [9] and references 
therein) . 



3. Unperturbed periodic counterpart: Matched quarter stack 

Without disorder (g a ^ = and r] a ,b = 0) all unit cells are identical since the parameters 
of a - and b -layers do not depend on the cell index n. Therefore, the array of bilayers is 
periodic with the period d which is the unperturbed size of a unit (a, b) cell, 

d = d a + d b . (3.1) 

As is known, the transmission through a periodic bilayer stack-structure is governed by 
the following dispersion relation [1], 

cos 7 = cos (p a cos ip b ( — " + — - J sin (p a sin </?&, (3.2) 



2 \Z b 

which defines the Bloch wave number k = j/d. We remind that ip a and ip b are the 
unperturbed phase shifts given by (12. lip . 

The dispersion relation (13. 2p specifies the band-structure in the dependence of the 
Bloch phase 7 on the wave frequency cu. Within the spectral bands where | C0S7I < 1, the 
solution j(cu) is real, therefore, electromagnetic waves can propagate through the bilayer 
stack. Otherwise, 7 is purely imaginary inside the spectral gaps where | COS7I > 1; here 
the waves are the evanescent Bloch states attenuated on the scale of the order of I7I -1 . 
As a result, inside spectral gaps the transmission is exponentially small for the samples 
with sufficiently large number N of unit (a, b) cells, N ^> I7I -1 . 

In what follows, we assume that without disorder the basic a- and b -layers are 
perfectly matched which means that their unperturbed impedances are equal, 

Z a = Z b . (3.3) 

In this case there are no waves reflected from the interfaces between matched layers. 
Therefore, in accordance with the dispersion relation (13.21) . such a stack-structure is 
equivalent to a homogeneous medium with the linear spectrum and mean refractive 
index n, 

1 1 _i _ n a d a + n b d b 

k = 7 jd = um/c, n= — — . (3.4) 

d a + d b 



Note that for this structure there are no gaps in the spectrum (I3.4p . and the Bloch phase 



is simply the total phase shift of the wave after passing any unit (a, b) cell, 

7 = <p a + ip b = u(n a d a + n b d b )/c. (3.5) 
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As we noted above, for a homogeneous RH-RH array all optic parameters are 
positive quantities. On the other hand, for a mixed RH-LH stack the parameters e a , 
fi a , and n a of a -layers are positive, whereas in b -layers £&, fib, and n& are negative. 
Therefore, the phase shift in the LH 6-layers is negative, ipb = —u\nb\db/c, however, due 
to their definitions both impedances, Z a and Z&, remain to be positive. As one can see, 
the only difference in the dispersion relation (13.21) is due to the sign before the second 
term. Specifically, for the RH-LH stack the sign "minus" in (13.21) has to be replaced 
by "plus", as compared with the RH-RH stack. Correspondingly, in the expressions 
(13. 4p and (13.51) for the RH-LH stack one has to make the substitution, rib ~^ ~ |^&|- 
As we show below, the change of the sign in the dispersion relation results in a strong 
change of properties of the localization length, therefore, in basic characteristics of the 
transmission. 

Let us now address the Hamiltonian map (12.1 2D - (12.171) . For a periodic (without 
disorder) bilayer stack the coefficients (12.141) are independent of the cell index n. 
Therefore, in line with (13. 3ft and (13.51) the Hamiltonian map (12.121) takes the following 
form 

Q n+ i = Q n cos 7 + P n sin 7, P n+i = -Q n sin 7 + P n cos 7. (3.6) 

Correspondingly, the map (12.171) in the radius-angle presentation transforms into the 
relations 

Rn+l = Rn, 6n+l = ®n ~ 7- (3-7) 

In spite of apparently simple forms of the unperturbed maps (13.61) . (13.71) their trajectories 
can be highly non-trivial. To show this, let us consider the matched quarter stack for 
which two basic a - and b -layers not only have the equal impedances (13.31) , but also the 
same lengths of optical paths, 

n a d a = \n b \d b . (3.8) 

For the homogeneous RH-RH array the refractive index of any b -layer is positive, nb > 0, 
therefore, the phase shifts are exactly the same, <p a = <pb, and the Bloch phase is always 
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non-zero, 7 = 2(p a . In this case the unperturbed map generates a circle in the phase 
space (Q, P) with a fixed radius R n . As for the angle 8 n , it changes by the Bloch phase 
7 in one step of the discrete time n, see Figure [TJ 

On the contrary, for the mixed RH-LH stack, the refractive index of b -layer is 
negative, nb < 0, therefore, <pb = —tp a - This means that the phase shift if a of the 
wave, gained in any RH layer, is canceled by the subsequent negative shift (fb in the 
next LH layer. Therefore, the phase shift 7 vanishes after passing every unit (a, b) cell, 
7 = 0. Together with the perfect transmission, this means that the structure consisting 
of a finite set of the RH-LH cells, has to be invisible for an observer. For such a kind 
of systems, the unperturbed trajectory in the phase space (Q, P) degenerates into a 
single stationary point, the fact which is drastically different as compared with the 
homogeneous RH-RH stack. 

Thus, when using the perturbation methods one has to take into account that in 
the zero order of perturbation the properties of trajectories of the Hamiltonian map 
remarkably depend on the type of a bilayer structure (fully conventional or having the 
metamaterial parts). Evidently, this fact should be reflected by the peculiarities of the 
Anderson localization. In particular, one can expect that the localization length for the 
RH-LH structure has a non-standard dependence on the disorder, as well as on control 
parameters of the underlying unperturbed system. 

4. Weak compositional disorder: Phase distribution 

Now we are in a position to turn on the disorder. As was shown in (9j [28], if the 
unperturbed impedances of two basic layers are equal, Z a = Zb, the localization 
length caused by the positional disorder, diverges. Indeed, the Anderson localization 
is originated from the effect of multiple wave reflections. However, the impedances 
are independent of the disorder in the layer thicknesses d an , db n , see definitions (12. lap . 
( 12.161) . (12.2 ap . ( 12.2 6p . Consequently, the matching condition ( 13. 3p remains valid even 
in the presence of positional disorder. Therefore, there is no any reflection from the 
interfaces even if these interfaces are randomly appear in the structure due to the 
positional disorder. Thus, for the structure with matched layers the positional disorder 
has an impact on transport characteristics only being accompanied by the compositional 
disorder, since the latter destroys the matching of layers. 

Taking into account this fact, as a first step it is reasonable to analyse the 
localization length contributed by the compositional disorder alone (g a ,b = 0). In this 
case all layers have constant thicknesses d a or db, so that the size d = d a + db of a unit 
(a, b) cell is also constant. The effect of positional disorder will be analyzed in the last 
part of the paper. In order to proceed, it is convenient to rewrite the definitions ( 13. 3p . 
( 13. 8 p for the unperturbed matched quarter stack in the form, 

Z a = Z b , ¥a = <f, <p h = ±y, (f = un a d a /c = u\n b \d b /c> . (4.1) 
Hereafter, the upper sign corresponds to the homogeneous RH-RH array, while the 
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lower sign is associated with the mixed RH-LH stack. In order to develop a proper 
perturbation theory, we assume the compositional disorder to be weak, 

a\ < 1 and K,V?) 2 < 1. (4.2) 

In doing so, we substitute (14. ip into (I2.14p with g ab = 0, and expand the coefficients 
A n , B n , C n , D n up to the second order in the compositional perturbation ^ 0) b(n) <C 1. 
Using the uncorrelated nature ( 12. 3 p of the disorder, from the exact perturbed #-map 
(12.17P one can obtain its quadratic approximation, 

9 n+ i ~ n = - 7 - Va(n)U(9 n ) T Vb(n)U(6 n - 7 /2) - o$W(9 n ). (4.3) 

Here the functions U(9) and W(9) are defined by 

U{9) = <p + sin ip cos(2# — <p) , 

W{9) = <^[cos(2# - 2<p) ± cos(2# - 2 7 )] 

+ sin ^[sin 9 sin(# - <p) ± sin(^ - 7 /2) sin(fl - tp - j/2)) ^ ' ' 

+ sin 2 tp sin(46 l — 2<p — 7) cos 7 . 

According to definition fl 2 . 1 8 [) of the localization length, one has to know the phase 
distribution that emerges as a result of iterations of the Hamiltonian map. When this 
distribution is uniform, the average in (12.181) can be performed relatively easy. However, 
in specific cases like the standard Anderson model with the energy at the band center or 
close to the band edges, the analytical derivation of the phase probability density p{6) 
is a highly non-trivial task. The same kind of problem occurs in our model with RH- 
LH stacks. In order to obtain the phase distribution, one has to derive the stationary 
Fokker-Plank equation for p(9). This can be done in the way described in [25 | |26 | |8| 127] . 
Specifically, we rewrite the map ( 14. 3 p in the continuum limit, 

d0 = -7dt - U(9)rj a (t)dt T U(9 - 7 /2)r/ b (t)dt - tfW(9)dt, (4.5) 

where random variables n a (t)dt, r)b(t)dt are replaced, respectively, with the Wiener 
variables dW a = r) a (t)dt and dWb = r] b (t)dt. As a result, expression ( 14. 5 p transforms to 
the stochastic Ito equation, 

d# = -U(6)dW a T U(9 - 7/2)dW 6 - [7 + rfW(d)] dt. (4.6) 

In accordance with ( 12. 3[) . the Wiener variables have the white- noise properties |25j : 

(dW a , 6 ) = 0, (dW a dW b ) = 2a 2 v 5 ab dt. (4.7) 

Following the theory of stochastic differential equations [25], one can readily associate 
the Ito equation ( 14. 6 p with the stationary Fokker-Plank equation for the probability 
density p(9), 



^ [t/ 2 (9) + U\9 - 7/2)] P (9)+2± 



1 + W(9) 



a 2 



p{9) = 0. (4. 



This equation should be complemented by the condition of periodicity and by the 
normalization condition, 

p(8 + ir)=p{9), fd0p(0) = l. (4.9) 
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Note that the integration in (14. 9 p is performed within the interval 7r due to the periodicity 
of p(9) with the period it. From the above equation (14. 8p one can see that its solution is 
sensitive to whether the Bloch phase 7 is non-zero (RH-RH array) or vanishes (RH-LH 
array) . 

5. RH-RH matched quarter stack 

In such a structure the unperturbed optic path lengths of two basic layers are equal, 
and the Bloch phase (13.51) is non-zero, 

n a d a = n b d b <p a = <p b = <p -> 7 = 2<p = 2um a d a /c. (5.1) 

In this case the term in the Fokker-Plank equation (14.81) containing 7/a 2 prevails over 
the others for any, even arbitrarily small, value of the phase shift tp. Indeed, under weak 
disorder conditions (14. 21) and small <p <C 1 one can get the relations a 2 W (9) / '7 ~ a 2 <C 1 
and a 2 U 2 (9) / '7 ~ a v {a r) p) <C 1, which allow to neglect all the terms with W and U in 
( 14.81) . As a result, the phase distribution p(9) is uniform within the first order of 
perturbation theory, 

p(9) = 1/tt. (5.2) 

The dimensionless inverse localization length d/Li oc is derived by differentiating 
the map ( 14 . 3 [) with respect to 9 and substituting the result into (I2.18p . with a further 
expansion of the logarithm within the quadratic approximation in disorder. Then, the 
subsequent averaging can be readily performed over random entries r] a ^(n) and over 
the angle 9 with the uniform distribution function (15. 2ft . After some algebra, one gets 

dj Li oc = A = a 2 sin 2 tp, tp = un a d a / c. (5.3) 

Note that in the #-map ( 14. 3 j) only linear terms in the compositional perturbations 
Va,b{ n ) contribute to the Lyapunov exponent, since the last quadratic term vanishes 
after averaging over 9. The final expression ( 15. 3ft is in a complete correspondence with 
the results previously obtained in the papers [291 [23J [21]. When the phase shift <p is 
small, the result (15. 3ft yields the asymptotics, 

d/L loc = \^o- 2 u 2 (n a d a /c) 2 for u<tic/n a d a . (5.4) 

This gives rise to the standard quadratic w-dependence A oc c 2 u 2 in the limit u — > 0. 

It should be noted that the localization length L[ oc (u) defined by (I5.3p . exhibits the 
Fabry-Perot resonances emerging when the phase shift p of the wave passing any layer 
is multiple to n, i.e., when the wave frequency u meets the condition 

oj / c = sir / n a d a s = 1,2,3,.... (5.5) 

At the resonances the factor sirup in expression (15.31) vanishes, thus resulting in the 
resonance divergence of the localization length Li oc and consequently, in the suppression 
of the localization. Remarkably, in contrast to the localization length Li oc , the Fabry- 
Perot resonances in the Lyapunov exponent A represent quite broad oscillations because 
of a smooth sin- function in expression ( 15.41) . 
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6. RH-LH matched quarter stack 

The principally different situation emerges for the mixed RH-LH matched quarter stack 
for which 

n a d a = -n b d b -> ip a = -<p b = <p ->■ 7 = 0. (6.1) 

In this case the Bloch phase vanishes independently of the value of the phase shift tp. As 
a result, the functions (14.41) turn out to be related to each other, W(9) = —U(8)U'(8), 
and the Fokker-Plank equation (14. 8 j) reduces to 

±U(0)±U(0)p(O)=O. (6.2) 

Solving this equation with the conditions (14. 9 p of periodicity and normalization, we get 
a highly nonuniform phase distribution, 

p(9) = -sin 2 tp/u(9). (6.3) 

Figure [2] displays a perfect agreement between analytical expressions (15. 2p . (16. 3p for the 
probability density p(9) and the corresponding numerical data obtained by the direct 
iteration of the exact map (12.121) . 

The Lyapunov exponent can be derived according to its definition (I2.18P and 
quadratic approximation (14.31) of the Hamiltonian #-map, taking into account the 
conditions (16.11) . In the second order of the perturbation theory in the weak 
compositional disorder r) a>b (n) one gets 

d/L loc = \ = 2a 2 r/ sin v?(cos(2# - <p)U(8)} = 0. (6.4) 

The averaging in (16.41) is performed over the angle 9 with the distribution function p{9) 
determined by (16. 3p . Since the denominator U(9) in (16.31) is the same as the coefficient 
in (16 .4p . we come to a very unexpected result: the Lyapunov exponent vanishes within 
the second order perturbation theory for any value of the phase shift ip [23] . This means 
that in order to derive a non- vanishing Lyapunov exponent, one has to go beyond the 
second order perturbation theory. For this, one has to obtain the expressions for both 
the #-map and the phase distribution p{9) in the next (fourth order) approximation, 
which is not a simple task. The problem is that the evaluation of high order terms in 
p{9) is not possible with the method [25], EHJ [HI [27] applied above. The reason is that the 
6*-map written within the fourth order approximation, cannot be reduced to the Fokker- 
Plank equation within the linear theory of differential stochastic equations. The only 
result which here can be drawn is that the Lyapunov exponent for the RH-LH matched 
quarter stack has to be proportional to er^, which is in contrast with the conventional 
quadratic dependence emerging in various disordered models. 

In order to proceed further, here we suggest another method that turns out to be 
very effective for deriving the expression for the Lyapunov exponent. The numerical 
data in Figure [2] give a hint of what can be done analytically. Specifically, for the RH- 
RH stack the trajectory in the phase space shown in Figure [2^, looks like it is created 
by the combination of two processes. The first one is a random increase of the radius 




Figure 2. (Color online) Phase trajectories generated by map (|2.12[) : (a) RH-RH 
matched quarter stack with unit-cell number N = 10 4 , cp = 27r/30, for zero disorder 
(solid circle), and for = 0.003 (scattered points); (b) RH-LH matched quarter stack 
with N = 10 6 , ip = 2tt/5, a* = 0.003. Phase distribution p(9): (c) for the RH-RH 
array from map (|2.12[) (histogram), and from expression (|5.2[) (horizontal line); (d) 
for the RH-LH array from map (|2.12[) (histogram), and from expression (|6.3I) (solid 
curve). Due to periodicity, p{9 + tt) = p{6), only the range ^ 9 < n is shown in (c) 
and (d) (after [54]). 



R n defined by (12.1 7ft . The second process is a uniform-like filling of a circle that is 
confirmed by the phase distribution p{9) generated by the exact Hamiltonian map. As 
one can see in Figure [2b, the phase distribution is uniform, the fact which has been 
already used for the analytical evaluation of the Lyapunov exponent. The right panel 
in Figure |2] demonstrates completely different effect: for the RH-LH stack the points fill 
an ellipse characterized by some aspect ratio, and rotated by some angle r with respect 
to the axes. This ellipse randomly increases in time, apparently keeping both the aspect 
ratio and the angle r. Thus, one can suppose that if to pass properly from the variables 
Q, P to new ones by both rotating the axes and rescaling ellipse into the circle, the 
trajectory in new variables will be of the same kind as those shown in Figure EK In 
such a way one can expect that the distribution of a new phase 9 n will be uniform, at 
least, approximately 

Following this idea, we make the linear transformation from "old" variables Q n , P n 
to the "new" ones Q n , P n , 

Q n = a 1/2 Q n cos r + a 1/2 P n sin r, 

~ (6-5) 
P n = -a~ 1/2 Q n sin r + a~ 1/2 P n cos r, 
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where the rotating angle r and rescaling a are the parameters to be specified. Note 
that in the new variables the expressions (12.12)) and (I2.16P - (I2.18P conserve their forms, 
however, with the new factors, 

A n = A n cos 2 r + (B n — C n ) sin r cos r + D n sin 2 r, 

B n oT x = B n cos 2 r — (A n — D n ) sin r cos r + C n sin 2 r, 

„ (6.6) 

C n a = C n cos 2 r + (A n — D n ) sin r cos r + £> n sin 2 r, 

D n = D n cos 2 r — (B n — C n ) sin r cos t + A n sin 2 r, 

which replace the old ones (12. 14ft . 

Now the distribution function p(9) for new phase 9 can be found starting from the 
quadratic expansion of the exact #-map (I2.17P with new coefficients (16. 6p . Taking into 
account the relations (16. ip . we come to the following equation, 

e n+ i - n = [f) a (n) - r] b (n)]V(9 n ) + a 2 V (9 n )V \9 n ) . (6.7) 

The stationary Fokker-Plank equation corresponding to this new #-map can be obtained 
in the same way as described before. Then, one gets 

iv(0)iv(0W0) = O. (6.8) 
d# dfl 

Note that the new Fokker-Plank equation (I6.8P differs from the old one (16.21) only in the 

function V(8) instead of U(6). The new function V(8) is defined by 

~ ~ ct ~ 

V(9) = siny9sin(2r — ip) sin 29 + — [ip — siny9cos(2r — tp)}(cos29 — 1) 

a- 1 

— [tp + sinv3cos(2r - <p)](cae29 + 1). (6.9) 

One can see that V(9) = —U{9) at a = 1 and r = 0, when new and old phases coincide. 

From equation (16. 8p with additional conditions similar to (I4.9p . one readily gets 
that the phase distribution is, indeed, uniform, p{9) = 1/ir, when the function V{9) is 
actually independent of the angle 9, i.e., when its derivative with respect to 9 vanishes, 

V'{9) = 0. (6.10) 

With this condition and definition (16. 9p . we can specify the values of r and a which 
determine the transformation from old to new phase-space variables. Also, one can 
obtain explicitly the ^-independent expression for the function V(9), 



r=|, a 2 = ^±^ ; vfy = V{) = Ji _ sin 2 (6n) 
2 ip — sin ip v 

The data presented in Figure [3] confirm the success of our approach: in the new variables 

the trajectory looks like a circle with the fluctuating radius. These data correspond to 

our expectation that the distribution of a new phase can be considered as uniform. 

The results (16. 6p . (16. lip allow us to calculate the Lyapunov exponent A for the 

mixed RH-LH matched quarter stack with the compositional disorder. For a weak 

disorder (I4.2p . the corresponding asymptotics of the #-map reads 

e n+1 -9 n = [ Va (n) - Vb (n)} V(ip) + \ [rjlin) - V 2 b (n)] Y(9 n ). (6.12) 
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Figure 3. (Color online) (a) Trajectory {Q ni P n ) generated numerically by the 
transformed map with factors (|6.6p and parameters (|6.11[) , for 7 = 0, ip = 2ir/5, 
ui = 0.02 and N = 10 7 ; (b) The corresponding distribution function p(9) (after [H]). 



Here the function Y{9) is defined as 

~ (sin -y? cosy?) sin ^ (2<p 2 - sin 2 cp) coscp - ipsimp ~ 

Y{d) = v&) + np] cos 29 ■ (6 - 13) 

Surprisingly, the fourth-order terms in the expansion of the #-map do not contribute in 
the fourth-order approximation for the Lyapunov exponent. This happens as a result of 
the averaging over new phase 9 with the uniform distribution. Such an average is in the 
spirit of the perturbation theory: the next fourth-order terms can be approximated by 
their average over the phase, with the phase distribution obtained in the second order 
approximation. 

Now we substitute the expression ( 16. 12ft into definition ( 12. 18ft for the Lyapunov 
exponent. Taking into account that A vanishes within the quadratic approximation in 
disorder, we expand the logarithm up to fourth-order terms in the perturbation. After, 
using the correlation properties ( I2.3P and averaging over 9 with uniform distribution 
function, we arrive at final result EE 



d 4 C + 2 [{2(f 2 - sin 2 ip) cos (p - up sin y] 2 

\ = a v — r^- . (6.14) 



Aoc 4 <y2 z - sin ip 



Here 



C = , 7 7 L ~ 3 (6-15) 

is the so-called excess kurtosis, the parameter which is specified by the form of the 
distribution of random variable r] a ,b{n). In general, ( can have any value within the 
range —2 ^ £ < 00. In particular, for the Gaussian and flat distributions it equals 
( = and £ = —6/5, respectively. 

Expression ( 16. 14p determines the asymptotics for small phase shift tp, 

Z^ = A ^ J^ + 2 K^ 8 for ^«l«^" 2 - (6-16) 
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Figure 4. (Color online) Localization length for RH-LH and RH-RH matched quarter 
stack of different lengths N versus the normalized wavelength 2n/tp = 2irc/n a d a uj, for 
af ? = 0.02. Full solid curve presents exact expression ()6.14|) . dashed-dotted line shows 
approximate dependence (|6.16|) . Lowest solid curve and dotted line correspond to the 
expressions (|5.3[) and (|5.4j) for RH-RH stack. 

As one can see, a quite surprising frequency dependence of the Lyapunov exponent, 
A oc cr^u 8 , emerges when uj — > 0. Thus, the dependence A oc w 6 numerically obtained 
for small values of u in Refs. [T4l H5] . can not be considered as the asymptotical result. 
From the analysis of the correct expression ( 16.1 4p it follows that the dependence A oc u> 6 
should be regarded as the intermediate one, apparently emerging due to not sufficiently 
large lengths iV over which the average of A was performed. 

Figure |4] demonstrates an excellent agreement between the numerical data for the 
localization length and analytical results derived above. The data have been obtained by 
the iteration of the exact map (I2.12p . with the use of definition (I2.18p . and compared with 
the analytical expressions ( I6.14p . ( 16 . 1 6[) for the RH-LH quarter stack with ( = —6/5, 
as well as f)5.3p . (15. 4p for the RH-RH structure. For sequences of unit-cell numbers 
N = 10 5 , 10 7 and 10 9 the data are obtained with the ensemble averaging over 100 
realizations of disorder, while for N = 10 12 only one realization is used. 

7. RH-LH matched quarter stack: Effect of positional disorder 

In the previous sections we presented a full analytical description of a nontrivial effect 
of compositional disorder in the mixed matched quarter stack with the LH material. On 
the other hand, as is noted in Section HJ the positional disorder itself in a bilayer stack 
with matched layers (Z a = Z^) does not result in the localization [91 [28] - Below we 
show that in the presence of an additional compositional disorder, one has to take into 
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account its non-trivial influence on the localization length. Specifically, we consider an 
interplay between both disorders, and show how the conventional frequency dependence 
of the localization length, L^ c oc u 2 , (for small u) emerges due to the influence of the 
positional disorder. 

Since we consider the mixed RH-LH matched quarter stack (14.11) . (16. ip . it is 
reasonable to treat from the beginning the transformed Hamiltonian map, which is 
described by equations (12121) . (12161) and Q2TTD with factors flEl} instead of (j2TT3j) . 
The angle r and rescaling parameter a are specified by expressions ( 16. lip . 

As before, in order to develop the perturbation theory we assume both 
compositional and positional disorders to be weak, compare with ( 14. 2ft . 

a 2 < 1, (or v (f) 2 < 1 and (a^) 2 < 1. (7.1) 

Note that the impedances (I2.2aj) . (12.2pp . depend only on the compositional disorder 
Va,b(n), whereas the phase shifts ( 12. 10 all . ( 12.10 6}) are randomized due to both 
compositional rj a>b {n) and positional g a ,b( n ) disorders. The latter is additionally 
incorporated into the layer thicknesses 1 12. Id) . ( 12 . 1 £>h . 

The expansion of the coefficients (16.61) and the recurrent relation (12. 1 7[) within the 
quadratic approximation in the disorders r] ajb (n), g a ,b( n ) yields the following ^-map, 

9 n+1 - 9 n = [rj a {n) - r) b (n)] V(<p) + [g a {n) - Q b (n)]G(9 n ) 

+\ [vl(n) - vl(n)} Y(9 n ) + \ [g 2 a {n) + o 2 b (n)] G(9 n )G'(9 n ). (7.2) 

This expression differs from (16 . 1 2[) in the second and fourth terms that take into account 
the effect of positional disorder. Functions V(<p) and Y{9) are respectively determined 
by (16. lip and (16.131) . As for the function G(9), it is given by 

~ ^(^-sin^cos2g) 

m = ' (7 - 3) 

The stationary Fokker-Plank equation corresponding to the #-map (I7.2p reads 

°1 V2 ^)&A°) + ^ 2 AG{9)tG{9)p{9) = 0. (7.4) 
a9 z a9 a9 

As one can see, in this equation the first term containing the variance a 2 , is responsible 

for the compositional disorder, and the second one (with a 2 ) arises due to positional 

disorder. As both disorders do not correlate with each other, the equations ( 17.21) and 

(I7.4p do not include the cross-correlated terms with the product a v a g . 

The solution to ( I7.4p satisfying the conditions of periodicity and normalization ( 14.91) 

is 

M = m [°IG\9) + a 2 V 2 (^)] " V2 , (7.5) 
with the normalization function /(y?) determined as 

r\<p) = ^ \a 2 e G 2 {9) + a 2 V 2 ^)} d9. (7.6) 
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Remarkably, when the positional disorder vanishes the distribution function p(9) reduces 
to the uniform one, 



p(9) ->■ 1/tt for a 2 Q ->■ 0. (7.7) 

This fact is in a complete accordance with the result obtained in the previous section. 

The dimensionless inverse localization length d/Li oc is obtained with the use of the 
map (17.21) in the same way as described above. Taking into account the white-noise 
properties (I2.3jl . the averaging over 9 with the probability density (17. 5ft . (I7.6jl yields 



d 



L 



A 



-o~ „ 



loc 



G(9)G"(9) ap\9) + a*V*(ip) 



-1/2 



d9 



C + 2 [{2ip 2 — sin 2 if) cos ip — ip sin ip]' 



4 



y? 2 — sin 2 if 



(7- 



From Figures and E] one can see that expression (17. 8p displays a very good 
agreement with the numerical data for the Lyapunov exponent obtained from the exact 
map ( 12.121) for the mixed RH-LH matched quarter stack. It is important to emphasize 
that this expression is applicable for any ratio between the variances <r 2 and <r 2 . Also, one 
should note that it is valid within a wide frequency range restricted to the requirement 
(17. ip of a weak disorder only. Remarkably, in the absence of compositional disorder, 
er 2 = 0, the Lyapunov exponent (17.81) vanishes even if the positional disorder is non- 
zero, a 2 7^ 0. This fact is a direct consequence of the layer matching (13.31) . and is in the 
accordance with the previously obtained results [9j 128] . 



u 




Figure 5. (Color online) Localization length versus the normalized wavelength 
2tt/(p = 2-Kc/n a d a ijj for RH-LH matched quarter stack randomized by both types of 
disorder with a\ = 0.01, o% = 0.00007 (left), and a 2 g = = 0.003 (right). The 
continuous curve presents numerical data for the structure length N = 10 7 with 
an ensemble averaging performed over 10 3 realizations of disorder. Large squares 
correspond to (|7.8[) and small squares stand as a reference to A oc ip 2 . 



The first term in (17. 8p prevails over the second one and correctly describes the 
localization length when the positional disorder predominates, <r 2 <C a 2 , or when both 
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Figure 6. (Color online) Same as in Figure ED for a 2 g = 0.00003 and a 2 = 0.02. The 
continuous curve presents numerical data, full squares correspond to expression (|7.8p . 
its first and second terms are shown by dashed-dotted and dashed curves, respectively. 

disorders are of the same order of magnitude, a 2 ~ a 2 , see Figure El The most important 
and interesting consequence of expression (17. 8p is that at any finite but sufficiently small 
strength a 2 of the compositional disorder, its second term is negligible. For this case we 
arrive at simple, however, non-trivial result, 

d/Li oc = A « a 2 sin 2 tp for a 2 < a\, (7.9) 

This asymptotics coincides with the Lyapunov exponent (15.31) for the RH-RH matched 
quarter stack and provides the conventional quadratic frequency dependence, A oc o~ 2 oj 2 
when u — > 0, see Figure El One should emphasize that the result ( 17. 9ft is valid when the 
positional disorder predominates. However, due to the layer matching ( 13. 3p . it does not 
contribute to the Lyapunov exponent. Thus, expression ( 17.91) turns out to be provided 
by small compositional disorder, instead of a large positional one. In spite of this fact, 
the presence of positional disorder plays the fundamental role, changing the abnormal 
octal frequency dependence (I6.16j) . L^ oc oc cr^u 8 , into the convetional quadratic one, 

Aoc 00 <7 ri Uj2 ^ wnen W — )■ 0. 

The second term in ( 17. 8ft coincides with the Lyapunov exponent ( 16.14p caused by 
the compositional disorder only. However, when the positional disorder is also present, 
this term can significantly contribute to the Lyapunov exponent if the compositional 
disorder prevails over the positional one, a 2 <C a 2 . At first glance, in this case one 
can neglect the terms with a 2 in the radicand of the first summand in ( 17. 8p and in 
the radicand of the normalization function ( 17.61) . In other words, one can suggest that 
the distribution function ( 17. 5p is uniform ( 17. 7p . However, a more careful analysis shows 
that this is not true at relatively low frequencies u, when the phase shift tp is sufficiently 
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small, y3 6 cx^ <C o 2 & <C ai. For such frequencies the Lyapunov exponent again obeys the 
standard quadratic cu-dependence originated from the first term in (17. 8p and is described 
by expression (17.91) . Nevertheless, when the condition a 2 g <C a 2 is met, the second 
term does contribute to the localization length in the regime of high and intermediate 
frequencies. This fact is clearly displayed in Figure El 

8. Conclusion 

In this paper we have studied two ID layered models which consist of both right-handed 
(RH) and left-handed (LH) materials for alternative a- and b -layers of thicknesses d a 
and db, respectively. The problem considered above is an analytical evaluation of the 
Lyapunov exponent whose inverse value is the localization length. As was found in a 
set of numerical studies [HI [15], for the RH-LH model a quite unexpected Anderson- 
type localization emerges, that has led to an intensive discussion in literature. In 
previous papers [23l [23] we have shown that the origin of this abnormal localization 
can be associated with a highly non-homogenous distribution of phase of the wave 
propagating along the bilayer structure. Unexpectedly, it was discovered [23] that with 
weak fluctuations of dielectric constants e an , Eb n the Lyapunov exponent vanishes in 
the second order of perturbation theory. This fact has shed light on the origin of the 
observed non-conventional localization. On the other hand, it was understood that the 
analytical treatment of this localization is of a very difficult task since one has to develop 
a proper perturbation theory beyond the second-order approximation. 

The problem of this unusual localization was rigorously solved with the use of 
a new method, and brief communication was published in [23]. Above, we have 
presented an analysis of the problem with many details that are important to understand 
our approach, together with a discussion concerning the mechanism of the observed 
phenomenon. 

We have to emphasize that originally the effect of abnormal localization was 
observed in [131 [15] for a quite specific bilayer model for which without disorder the RH- 
LH array represents an ideal mixed stack. Specifically, it was assume that e a — \i a — 1, 
£b = f'b = — 1) hence the impedances are trivially equal, Z a — Z b — 1. The additional 
assumption was that the layer thicknesses are also equal, d a = db- In the present 
study we consider a more general model of the so-called matched quarter stack for 
which the condition of equal optical pass lengths, n a d a = \n b \db, and the equality of the 
impedances, Z a = Zb, are independent of each other. Thus, our results can be applied 
to a more general class of physical systems. 

In our approach we use the reduction of the underlying problem of a quantum 
localization to the study of the classical Hamiltonian map of a quite transparent 
structure. In this way the evaluation of the localization length can be performed in 
terms of local instability of trajectories in the phase space of the corresponding classical 
model. Then, with the transformation to radius-angle variables in the phase space, the 
equations allowing to evaluate the Lyapunov exponent take a quite simple form, and the 
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problem is reduced to the analytical treatment of the angle distribution which emerges 
after a sufficiently long time of iteration of the classical map. It should be noted that 
the angle in this classical map is nothing but the phase of wave function in the quantum 
problem. Therefore, even with the reduction of the quantum model to the classical one 
the physical correspondence between two representations remains transparent. 

Our approach allows one to resolve the problem of the Anderson localization for 
bilayer RH-LH arrays with the matched layers, and derive the localization length. As we 
show, the peculiarity of this model is entirely due to the zero value of the unperturbed 
Bloch phase that happens when the unperturbed impedances are equal. We rigorously 
derive the expression for the Lyapunov exponent which appears to be defined by the 
fourth order of perturbation theory, and which is valid for any value of frequency u. Our 
results prove that for small u the localization length is enormously large, L\ oc oc a~ 4 u~ s . 
The generalization of our result is due to the ^-map (14.31) . according to which the 
same dependence for L ioc is expected to occur when the effective unperturbed refractive 
index n = (n a d a — |n&|<4)/(d a + c4), see (13. 4p . is sufficiently small, n <C a 2 . Thus, the 
obtained results show that the non-conventional localization observed numerically in 
Refs. [TH [151 EEl EH] is ver y fragile with respect to the choice of model parameters, and 
can be hardly observed experimentally. 

It is interesting that the discussed effect of anomalous localization can be considered 
as a kind of strongly correlated disorder, one part of which is embedded into the 
randomized phase shifts ip an , ipb n and the other is absorbed by the impedances Z an , 
Zb n . As one can see, the compositional disorder, even if it is of a white-noise type, 
results in strong correlations between fluctuations of phase shifts and impedances. To 
the best of our knowledge, this phenomenon of a natural emergence of correlated disorder 
from the white-noise disorder was never discussed in literature. 

A special attention in our study was paid to the model with two kinds of disorder. 
Specifically, we have considered the case when in addition to the compositional disorder 
with the fluctuating values of dielectric permittivities e an , Eb n , the positional disorder 
is also included (therefore, the thicknesses d an and/or d bn also slightly fluctuate). It is 
extremely important that the positional disorder randomly affects the phase shifts only 
and do not contribute to the impedances. As one can see, the positional disorder plays 
a fundamental role for the destroy of strong correlations between the disordered phase 
shifts and impedances, therefore, for the recovering of the conventional dependence 
Lf oc oc cr^u 2 when oj — > 0. 

The positional disorder in both thicknesses d an and db n has been earlier analysed 
in [28]. Specifically, the general case of correlated disorder was analytically studied 
by assuming any kind of statistical correlations in thicknesses of a- and & -layers, as 
well as the inter-correlations between the fluctuations of two thicknesses. It was shown 
that for any ratio between d a and d^ the localization length is governed by the unique 
expression, no matter whether the structure represents RH-RH or RH-LH array. The 
main result is that for such structures the phase distribution is flat, therefore, the 
analytical expressions for the Lyapunov exponent can be obtained within the standard 
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perturbation theory. The same statement stems from the analysis of RH-RH and RH- 
LH arrays with fluctuating refractive indices [29], however, only when n a d a 7^ |n&|cZ& in 
the latter case. These results also indicate that the abnormal localization emerges in a 
very specific situation, and not only due to the inclusion of left-handed materials into 
the structure. 

Finally, we would like to note that the effect of vanishing of the unperturbed Bloch 
phase shift 7 in the #-map (I4.3p . is a typical effect occurring at band edges in the 
tight-binding Anderson model as well as in the Kronig-Penney models. Indeed, in 
these models after one period of perturbation the unperturbed Bloch phase 7 vanishes 
when approaching the band edges. This results in a very specific non- homogeneous 
distribution of the wave phase, therefore, to a non-standard expression for the Lyapunov 
exponent. However, in the considered model of the mixed RH-LH matched quarter 
stack, the peculiarity is that vanishing Block phase emerges independently of the value 
of frequency u. This is a principal difference as compared with the localization occurring 
in the Anderson and Kronig-Penney models for which the zero Bloch phase emerges at 
band edges, therefore, only for specific values of frequency. 

Acknowledgments 

We gratefully acknowledge the support of the SEP-CONACYT (Mexico) grant CB-2011- 
01-166382, as well as the support from VIEP grant EXC08-G of the BUAP (Mexico). 

References 

[1] Markos P and Soukoulis C M 2008 Wave Propagation. From Electrons to Photonic Crystals and 

Left-Handed Materials (Princeton, NJ: Princeton University Press) 
[2] McGurn A R, Christensen K T, Mueller F M and Maradudin A A 1993 Phys. Rev. B 47 13120 

Nau D et al 2007 Phys. Rev. Lett. 98 133902 

Ponomarev IV et al 2007 Phys. Rev. B 75 205434 
[3] Baluni V and Willcmsen J 1985 Phys. Rev. A 31 3358 

Vinogradov A P and Merzlikin A M 2003 Physica B 338 126 

Nascimento E M, dc Moura F A B F and Lyra M L 2007 Phys. Rev. B 76 115120 

Luna-Acosta G A et al 2009 Phys. Rev. B 80 115112 
[4] Bcllani V et al 1999 Phys. Rev. Lett. 82 2159 
[5] Czycholl G, Kramer B and MacKinnon A 1981 Z. Phys. B 43 5 

[6] Thoulcss D J, in Ill-condensed Matter, cd. G Toulouse, R Balian, Amsterdam, North-Holland, 
1979. 

[7] Kappus M and Wegncr F 1981 Z. Phys. B 45 15 

Dcrrida B and Gardner E 1984 J. Physique 45 1283 
[8] Izrailcv F M, Ruffo S and Tessicri L 1998 J. Phys. A: Math. Gen. 31 5263 
[9] Izrailcv F M, Krokhin A A and Makarov N M 2012 Phys. Rep. 512 125 
[10] Kuhl U, Izrailcv F M, Krokhin A A, and Stockmann H-J Appl. Phys. Lett. 2000 77 633 

Krokhin A A, Izrailcv F M, Kuhl U, Stockmann H-J, and Ulloa S E, 2002 Physica E 13 695 
Kuhl U, Izrailev F M, and Krokhin A A, 2008 Phys. Rev. Lett. 100 126402 
[11] Dictz O et al 2012 Phys. Rev. B 86 201106(R) 
[12] Shelby R A, Smith D R and Schultz S 2001 Science 292 77 



Non-conventional Anderson localization in a matched quarter stack 



23 



[13] Gorkunov M V, Gredeskul S A, Shadrivov I V and Kivshar Yu S 2006 Phys. Rev. E 73 056605 

Dong Y and Zhang Z 2006 Phys. Lett. A 359 542 
[14] Asatryan A A et al 2007 Phys. Rev. Lett. 99 193902 
[15] Asatryan A A et al 2010 Phys. Rev. B 81 075124 
[16] Asatryan A A et al 2010 Phys. Rev. B 82 205124 
[17] Mogilevtscv D et al 2010 Phys. Rev. B 82 081105(R) 

Mogilevtsev D et al 2011 Phys. Rev. B 84 094204 
[18] Maurel A, Ourir A, Mcrcier J-F and Pagneux V 2012 Phys. Rev. B 85 205138 
[19] del Barco O and Ortmio M 2012 Phys. Rev. A 86 023846 
[20] Mogilevtscv D et al. 2010 Phys. Rev. B 81 047601 

Asatryan A A et al 2012 Phys. Rev. B 85 045122 
[21] Reyes-Gomez, Bruno- Alfonso A, Cavalcanti S B, and Olivcira L E 2012 Phys. Rev. B 85 195110 
[22] Gredeskul A et al 2012 Low Temp. Phys. 38 570 [in Russian: 2012 Fiz. Nizk. Temp. 38 728] 
[23] Torrcs-Hcrrera E J, Izrailev F M and Makarov N M 2011 Low Temp. Phys. 37 957 [in Russian: 

2011 Fiz. Nizk. Temp. 37 1201] 
[24] Torres-Herrera E J, Izrailev F M and Makarov N M 2012 EPL 98 27003 
[25] Gardiner C W 2004 Handbook of Stochastic Methods 3rd Edition (Berlin: Springer) 
[26] Izrailev F M, Kottos T and Tsironis G 1995 Phys. Rev. B 52 3274 
[27] Hernandez-Herrejon J C, Izrailev F M and Tessieri L 2008 Physica E 40 3137 
[28] Izrailev F M and Makarov N M 2009 Phys. Rev. Lett. 102 203901 

[29] Izrailev F M, Makarov N M and Torres-Herrera E J 2010 Physica B: Condensed Matter 405 3022 



